Assignment 1 STAT5900F

Colon dataset

Before you start. This document has many graphic visualization and statistical test done. The aim is to obtain a better view of the data and think above this first task.

Some visualizations show the same result but different options of view are chosen to understand better the distribution and relation of the data.

Source of data: The data are described in Alon et al. (1999) and can be freely downloaded from http://microarray. princeton.edu/oncology/affydata/index.html.

first we load some libraries. And try to better undersand our dataset. (dimenssion and type of variables)

# load plsgenomics library
library(plsgenomics)

For any news related to the 'plsgenomics' package (update, corrected bugs), please check http://thoth.inrialpes.fr/people/gdurif/
C++ based sparse PLS routines will soon be available on the CRAN in the new 'fastPLS' package.

# load data set
data(Colon)
# how many samples and how many genes ?
dim(Colon$X) # 2000 genes from 62 samples

[1]   62 2000

# how many samples of class 1 and 2 respectively ?
sum(Colon$Y==1) # we have 22 samples of class 1

[1] 22

sum(Colon$Y==2) # we have 40 samples of class 2

[1] 40

class(Colon)
class(Colon$X)
typeof(Colon$X)
dim(Colon$X)

class(Colon$Y)
typeof(Colon$Y)
table(Colon$Y)

length(Colon$Y)
class(Colon$gene.name)
typeof(Colon$gene.name)

expression.mat<-t(Colon$X)
#head(expression.mat)# expression matrix has 62 genes

# visualize 6 histograms of the first 6 genes
par(mfrow=c(2,3))# prepare a window to visualize 6 histograms for the first 9 genes
for (i in 1:6)
{
hist(expression.mat[i,],freq=FALSE,main= paste("Expression values of Gene", i),xlab="Expression values",ylab="Relative Frequency")
lines(density(expression.mat[,i]))
}

par(mfrow=c(1,1))# return in one plot window

# Boxplot
boxplot(t(expression.mat[1:10,]),xlab="Genes 1-10",ylab="Expression values")

boxplot(t(expression.mat[11:20,]),xlab="Genes 11-20",ylab="Expression values")

boxplot(t(expression.mat[1:70,]),xlab="Genes 1-70",ylab="Expression values")# boxplot of 70 genes

boxplot(t(expression.mat),xlab="Genes",ylab="Expression values")# boxplot of all genes

# Correlations
xnew<-t(expression.mat[1:9,])
colnames(xnew)<-c("Gene 1","Gene 2","Gene 3", "Gene 4", "Gene 5", "Gene 6", "Gene 7","Gene 8","Gene 9")
pairs(xnew,pch=16,col="Blue")# we observe that Gene 2 and Gene 3 are correlated

par(mfrow=c(3,3))# prepare a window to visualize QQplot for the first 9 genes

for (i in 1:9)
{
qqplot(y=xnew[,i],x=rnorm(1000),pch=16,col="blue",ylab="Sample quatiles",xlab="Theoretical quantiles of standard normal",main=paste("Q-Q plot of Gene", i))
qqline(y=xnew[,i])
}

heatmap(expression.mat)

heatmap(expression.mat,Colv=NA)

heatmap(t(expression.mat),main="Original data")

NA
NA

Reference: https://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html

Before starting we have to install these packages limma and ComplexHeatmap,biobase from Bioconductor:

# if (!require("BiocManager", quietly = TRUE))
 # install.packages("BiocManager")

 #BiocManager::install("ComplexHeatmap")
 #BiocManager::install("limma")
 #BiocManager::install("Biobase")

library(ComplexHeatmap)

Loading required package: grid
========================================
ComplexHeatmap version 2.12.1
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite either one:
- Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
    genomic data. Bioinformatics 2016.
- Gu, Z. Complex Heatmap Visualization. iMeta 2022.


The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================

Heatmap(as.matrix(t(expression.mat)),row_split = factor(Colon$Y),cluster_columns = ,border=T,heatmap_legend_param = list(title = ""), cluster_row_slices = F)

The automatically generated colors map from the 1^st and 99^th of the
values in the matrix. There are outliers in the matrix whose patterns
might be hidden by this color mapping. You can manually set the color
to `col` argument.

Use `suppressMessages()` to turn off this message.

library(limma)
library(Biobase)
my_set = ExpressionSet(assayData=expression.mat,phenodata=Colon$Y)
plotDensities(my_set,legend="topright")

my_set = ExpressionSet(assayData=expression.mat[1:20,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright",main="First 20 genes")

my_set = ExpressionSet(assayData=expression.mat[1:50,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 50 genes")

my_set = ExpressionSet(assayData=expression.mat[1:100,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 100 genes")

my_set = ExpressionSet(assayData=expression.mat[1:500,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright",main="First 500 genes")

my_set = ExpressionSet(assayData=expression.mat[1:1000,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 10^3 genes")

my_set = ExpressionSet(assayData=expression.mat[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 2000 genes ")

Estimated variance

est_var <- apply(Colon$X, 2, var)


hist(est_var, prob = TRUE, main = "Estimated variances of all genes", xlab = "Estimated Variances", freq = FALSE, breaks = 50)
lines(density(est_var), lwd = 2)

NA
NA
NA

Normalized gene expression matrix

Approach 1

Reference: https://www.frontiersin.org/articles/10.3389/fgene.2019.00400/full After applying the function it is not seen too much change in the distribution of the data.

# Install R package
#install.packages("NormExpression");
#install.packages("dendextend");

# Load libraries 
library(NormExpression);
library(dendextend);

---------------------
Welcome to dendextend version 1.16.0
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
You may ask questions at stackoverflow, use the r and dendextend tags: 
     https://stackoverflow.com/questions/tagged/dendextend

    To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------


Attaching package: ‘dendextend’

The following object is masked from ‘package:stats’:

    cutree

# Using the parameters to produce the TU normalization factor
tu <- getFactors(data = t(Colon$X), method = "TU", pre_ratio=1, lower_trim=0.2, upper_trim=0.6);

Parameters range 0.2 - 0.6 ...identified too few ubiquitous genes ( 0 ), trying range 5-95  instead 

TU_col.matrix <- getNormMatrix(data =t(Colon$X), norm.factors = tu);
head(TU_col.matrix)

              1        2         3         4        5        6        7        8
Gene 1 9839.864 7775.411  6581.322  8048.674 4557.638 4747.677 8768.282 5195.505
Gene 2 6264.307 5701.185 11991.043 10080.780 5212.463 3708.107 4649.614 4070.246
Gene 3 4884.074 4143.362  9237.904  7674.213 4798.070 2962.309 4559.165 3701.538
Gene 4 4656.709 3154.672  8095.083  5122.594 4886.738 5811.494 8113.627 5697.014
Gene 5 2288.746 1709.815  2006.810  2580.407 3077.744 3423.568 3027.799 2391.092
Gene 6 6051.326 4725.788  2704.585  2745.247 4123.724 3165.145 4645.448 3647.437
              9       10       11        12       13       14       15        16
Gene 1 7296.218 4605.567 6916.785 10198.930 5375.243 7748.276 4707.026  8870.686
Gene 2 5487.193 4141.867 1488.725  6053.175 6214.729 7360.607 4995.598  5705.758
Gene 3 3693.809 2899.286 1082.230  4348.428 5741.591 6333.759 4108.151  4108.834
Gene 4 7422.000 5934.271 5230.074  6308.047 4005.629 4422.227 3078.209  4011.023
Gene 5 4180.782 2300.276 2740.080  2526.576 3817.851 4621.213 3014.960  9406.395
Gene 6 1217.623 4760.809 2473.468  3510.340 4382.479 6558.102 4272.751 12195.613
              17        18       19       20       21       22       23       24
Gene 1 14605.810 10144.130 7819.326 6205.795 6662.269 5568.276 7242.968 8583.538
Gene 2  4234.713  6750.863 5590.864 2464.604 5798.788 5299.788 3803.477 9500.019
Gene 3  3310.400  4881.235 4899.417 1800.673 4506.990 4231.798 2769.266 7152.250
Gene 4  5266.622  4947.110 6456.888 7313.869 4255.826 3221.781 3031.692 5316.998
Gene 5  6953.209  6376.588 3580.912 3331.894 1387.575 1318.164 7517.709 9371.252
Gene 6  5959.528  2808.811 1895.208 2880.278 6272.745 4964.393 7768.621 4876.405
             25       26       27       28       29       30       31       32
Gene 1 7473.582 7528.110 7655.893 6605.050 7083.159 7582.618 5125.371 6074.719
Gene 2 4591.978 4750.331 4929.246 6617.728 5052.686 4420.438 5257.057 1676.830
Gene 3 4119.789 4476.149 4523.814 5819.247 4114.954 3350.963 4554.237 1538.852
Gene 4 3188.727 3452.345 3291.601 3341.562 4243.503 4038.357 3631.810 3840.575
Gene 5 3310.470 1974.493 2037.854 2422.091 1352.836 1808.892 2501.083 2357.996
Gene 6 5175.153 4546.428 5465.097 4913.106 5402.603 5235.904 4830.048 3304.968
             33       34       35       36       37       38       39       40
Gene 1 8470.879 6043.356 6316.759 5595.525 3992.116 5343.519 4363.859 8653.000
Gene 2 2773.544 7132.374 4140.535 3856.801 7413.494 3084.937 6065.678 5216.300
Gene 3 2437.673 5920.998 3623.296 3634.545 5667.255 2834.632 4895.064 4351.217
Gene 4 4219.503 3234.159 3609.084 2037.313 5300.192 1912.300 1629.772 3066.051
Gene 5 3531.888 2452.998 1834.420 3927.101 1073.212 4838.643 2363.131 3771.462
Gene 6 7222.493 4077.580 4579.722 2801.599 4159.132 4567.388 5828.242 7277.141
             41       42       43       44       45       46       47       48
Gene 1 5921.792 4678.451 5807.956 7207.930 4093.071 6662.282 3119.976 4255.698
Gene 2 8831.477 3929.019 1790.503 1662.623 2891.615 1693.355 2498.088 2786.928
Gene 3 8038.527 3081.307 1481.889 1500.682 2713.414 1531.501 2011.111 2454.858
Gene 4 2119.561 2984.300 2762.958 2658.916 2126.183 2660.821 1646.150 1758.747
Gene 5 2918.967 1799.268 1369.471 2330.601 1360.582 2372.147 1920.372 1212.969
Gene 6 3457.447 1472.184 3367.695 2543.795 4735.532 2267.760 3397.162 1847.481
             49       50       51       52       53       54       55       56
Gene 1 7048.964 5944.613 6338.834 6331.638 7712.292 3288.024 5979.722 5066.804
Gene 2 4143.003 4463.275 4384.165 4516.564 6447.877 3923.083 3506.673 3259.570
Gene 3 2955.314 3551.252 3784.988 3991.300 6528.478 3609.653 2498.275 2378.372
Gene 4 2764.934 2838.009 2935.844 1965.263 2465.702 2491.357 2942.947 2920.153
Gene 5 1448.224 3353.343 2946.598 4415.614 4201.060 2401.942 1518.474 2325.745
Gene 6 4537.555 5976.121 6702.973 3585.764 4682.331 2679.671 1996.398 3232.997
             57       58       59       60       61       62
Gene 1 4808.376 8755.587 7638.355 6125.205 7260.405 7206.305
Gene 2 6459.984 7349.936 5720.557 3159.807 4664.292 3524.000
Gene 3 5112.060 6459.995 5015.061 2329.239 3602.677 2631.201
Gene 4 3889.023 2760.340 2637.977 2388.598 3706.799 3370.217
Gene 5 2828.418 2765.420 4063.208 1452.638 2090.693 2319.156
Gene 6 2885.635 6104.776 9407.301 3979.457 4589.234 5585.658

my_set = ExpressionSet(assayData=TU_col.matrix[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 2000 genes getNormMatrix ")

Approach 2

Using logarithm of base 2

log2_col.matrix <- log2(expression.mat)
head(log2_col.matrix)

         1        2        3        4        5        6        7        8
1 13.06834 13.16180 11.90151 12.60882 11.65747 11.29366 12.79900 11.97610
2 12.41686 12.71414 12.76702 12.93360 11.85114 10.93712 11.88381 11.62395
3 12.05779 12.25368 12.39070 12.54009 11.73163 10.61316 11.85547 11.48696
4 11.98902 11.86037 12.20018 11.95694 11.75805 11.58535 12.68705 12.10904
5 10.96426 10.97672 10.18804 10.96767 11.09105 10.82194 11.26497 10.85651
6 12.36696 12.44344 10.61854 11.05701 11.51313 10.70871 11.88252 11.46572
         9       10       11       12       13       14       15       16
1 13.18777 12.36400 13.86074 13.79086 12.12577 12.28344 12.26262 12.45821
2 12.77669 12.21091 11.64471 13.03821 12.33513 12.20939 12.34846 11.82159
3 12.20573 11.69632 11.18464 12.56101 12.22089 11.99263 12.06630 11.34790
4 13.21243 12.72969 13.45747 13.09771 11.70146 11.47434 11.64990 11.31314
5 12.38439 11.36243 12.52485 11.77770 11.63220 11.53784 11.61995 12.54281
6 10.60469 12.41183 12.37717 12.25213 11.83118 12.04284 12.12297 12.91746
        17       18       19       20       21       22       23       24
1 12.80260 12.24829 12.39403 12.10737 12.86008 12.77219 12.04128 10.90289
2 11.01639 11.66079 11.91006 10.77511 12.65982 12.70090 11.11202 11.04924
3 10.66113 11.19296 11.71960 10.32229 12.29623 12.37623 10.65421 10.63971
4 11.33101 11.21230 12.11783 12.34439 12.21351 11.98282 10.78482 10.21193
5 11.73181 11.57850 11.26731 11.21010 10.59663 10.69349 12.09499 11.02956
6 11.50933 10.39568 10.34934 10.99996 12.77317 12.60658 12.14236 10.08713
        25       26       27       28       29       30       31       32
1 13.11398 12.53501 12.50720 13.25370 13.70182 13.66375 13.09560 12.29118
2 12.41130 11.87075 11.87200 13.25647 13.21448 12.88524 13.13220 10.43410
3 12.25475 11.78498 11.74817 13.07097 12.91830 12.48563 12.92515 10.31021
4 11.88516 11.41029 11.28942 12.27066 12.96268 12.75482 12.59863 11.62968
5 11.93922 10.60420 10.59768 11.80639 11.31341 11.59616 12.06050 10.92592
6 12.58378 11.80745 12.02088 12.82677 13.31108 13.12949 13.00998 11.41300
        33       34       35       36       37       38       39       40
1 12.75338 13.02721 12.31643 12.15310 12.41652 12.03669 11.83636 13.15610
2 11.14260 13.26624 11.70707 11.61623 13.30952 11.24414 12.31142 12.42593
3 10.95638 12.99770 11.51455 11.53060 12.92201 11.12206 12.00208 12.16432
4 11.74795 12.12525 11.50888 10.69550 12.82541 10.55421 10.41542 11.65928
5 11.49131 11.72641 10.53257 11.64229 10.52130 11.89350 10.95145 11.95803
6 12.52337 12.45958 11.85251 11.15508 12.47564 11.81027 12.25382 12.90627
        41       42       43       44       45       46       47       48
1 11.89144 12.59671 13.48276 13.83782 13.20516 13.81060 12.76308 12.62177
2 12.46806 12.34485 11.78510 11.72169 12.70386 11.83447 12.44237 12.01105
3 12.33234 11.99422 11.51217 11.57385 12.61209 11.68953 12.12954 11.82802
4 10.40917 11.94807 12.41094 12.39907 12.26025 12.48645 11.84064 11.34692
5 10.87086 11.21809 11.39835 12.20894 11.61620 12.32078 12.06293 10.81092
6 11.11511 10.92864 12.69649 12.33522 13.41551 12.25585 12.88587 11.41794
        49       50       51       52       53       54       55       56
1 12.63434 12.79790 12.74616 13.50258 12.90439 12.02835 12.49615 12.14434
2 11.86761 12.38442 12.21425 13.01523 12.64605 12.28312 11.72617 11.50795
3 11.38024 12.05464 12.00223 12.83686 12.66397 12.16299 11.23700 11.05324
4 11.28417 11.73120 11.63572 11.81472 11.25922 11.62806 11.47333 11.34931
5 10.35121 11.97192 11.64099 12.98261 12.02798 11.57533 10.51869 11.02096
6 11.99884 12.80553 12.82675 12.68228 12.18445 11.73318 10.91347 11.49614
        57       58       59       60       61       62
1 12.18402 12.27966 13.15361 12.71652 12.60609 12.86728
2 12.61000 12.02719 12.73650 11.76160 11.96769 11.83523
3 12.27237 11.84099 12.54662 11.32163 11.59511 11.41374
4 11.87788 10.61430 11.61978 11.35793 11.63621 11.77086
5 11.41846 10.61695 12.24297 10.64044 10.81002 11.23162
6 11.44736 11.75939 13.45413 12.09434 11.94429 12.49975

my_set = ExpressionSet(assayData=log2_col.matrix[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 2000 genes->log2 ")

boxplot(t(log2_col.matrix[1:70,]),xlab="Genes 1-70",ylab="Expression values", main="Log2 transformed data")# boxplot of 70 genes

boxplot(t(log2_col.matrix),xlab="Genes",ylab="Expression values",main="Log2 transformed data")# boxplot of all genes

heatmap(t(log2_col.matrix),Colv=NA)

heatmap(t(log2_col.matrix))

library(ComplexHeatmap)
Heatmap(as.matrix(t(log2_col.matrix)),row_split = factor(Colon$Y),cluster_columns = ,border=T,heatmap_legend_param = list(title = ""), cluster_row_slices = F)

#Approach 3 Z scores

Using Z scores

score_col.matrix <-(expression.mat-mean(expression.mat))/sd(expression.mat) # by mean and sd of expression matrix
head(score_col.matrix[,1:5])

          1         2        3         4        5
1 11.628906 12.445549 4.861345  8.300372 4.015524
2  7.194812  8.972454 9.328798 10.540857 4.674878
3  5.483167  6.364028 7.055202  7.887512 4.257619
4  5.201208  4.708560 6.111438  5.074243 4.346901
5  2.264666  2.289284 1.083625  2.271372 2.525394
6  6.930691  7.339244 1.659860  2.453115 3.578609

my_set = ExpressionSet(assayData=score_col.matrix [1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend="topright", main="First 2000 genes-> z score ")

boxplot(t(score_col.matrix[1:62,]),xlab="Genes 1-62",ylab="Expression values", main="Z score transformed data")# boxplot of 62 genes

boxplot(t(score_col.matrix),xlab="Genes",ylab="Expression values",main="Z score transformed data")# boxplot of all genes

heatmap(score_col.matrix,Colv=NA)

heatmap(t(score_col.matrix))

NA
NA

library(ComplexHeatmap)

Loading required package: grid
========================================
ComplexHeatmap version 2.12.1
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite either one:
- Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
    genomic data. Bioinformatics 2016.
- Gu, Z. Complex Heatmap Visualization. iMeta 2022.


The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================

Heatmap(as.matrix(t(score_col.matrix)),row_split = factor(Colon$Y),cluster_columns = ,border=T,heatmap_legend_param = list(title = "Z-score"), cluster_row_slices = F)

The automatically generated colors map from the minus and plus 99^th of
the absolute values in the matrix. There are outliers in the matrix
whose patterns might be hidden by this color mapping. You can manually
set the color to `col` argument.

Use `suppressMessages()` to turn off this message.

my_set = ExpressionSet(assayData=expression.mat[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend=FALSE, main="2000 genes original")

my_set = ExpressionSet(assayData=TU_col.matrix[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend=FALSE, main="2000 genes getNormMatrix ")

 my_set = ExpressionSet(assayData=log2_col.matrix[1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend=FALSE, main="2000 genes->log2 ")

my_set = ExpressionSet(assayData=score_col.matrix [1:2000,],phenodata=Colon$Y)
plotDensities(my_set,legend=FALSE, main="2000 genes-> z score ")

End 3 approach of transforming

Extracting the table of differential genes

So, we decide to go with log2 transformation data and work with the transformed data. We can now have a look at the result table that contains all information on all genes (p-value, fold changes, etc). ## False discovery rates

# Recall our notations:
expression.mat<-t(Colon$X)
log2_col.matrix <- log2(expression.mat)
# as.matrix(t(log2_col.matrix))# log2 matrix transformed

library(limma)
library(Biobase)
my_set <- ExpressionSet(assayData = log2_col.matrix, phenodata = Colon$Y)

design <- model.matrix(~ as.factor(Colon$Y))
colnames(design) <- c("Intercept", "Tumor-Normal")
mod1 <- lmFit(my_set, design = design)
head(mod1$sigma^2)

   Gene 1    Gene 2    Gene 3    Gene 4    Gene 5    Gene 6 
0.3722873 0.4389969 0.4460575 0.4828473 0.3829682 0.6324110 

mod2 <- eBayes(mod1)
head(mod2$s2.post)

   Gene 1    Gene 2    Gene 3    Gene 4    Gene 5    Gene 6 
0.5092709 0.5588206 0.5640650 0.5913912 0.5172043 0.7024821 

mod2$s2.prior

[1] 0.904814

# topTable(mod2, )

results <- decideTests(mod2[,"Tumor-Normal"]) # 
final <- topTable(mod2, coef=2, adjust="BH",sort="P",n=Inf) # sorted by p-value
head(final,n=10)

summary(results)

       Tumor-Normal
Down             29
NotSig         1857
Up              114

vennDiagram(results) #how many significant

log2(1857/143)

[1] 3.698887

library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:Biobase’:

    combine

The following objects are masked from ‘package:BiocGenerics’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

# distribution of not adjusted p-values
hist(final$P.Value, col="lightblue", main = "NOT adjusted p-value distribution")

# distribution of adjusted p-values
hist(final$adj.P.Val, col="lightblue", main = "Adjusted p-value distribution")

Differential Expression Analysis by using two-sample t-test

p_twosample<-NULL
t_twosample<-NULL

for (i in 1:ncol(log2(Colon$X))){
temp<-t.test(log2(Colon$X)[,i]~Colon$Y , alternative = "two.sided", var.equal = TRUE)  
t_twosample[i]<-temp$statistic
p_twosample[i]<-temp$p.value
}

par(mfrow=c(1,2))
hist(final[,3],main="Moderated t-statistics",breaks=15,xlab="Moderated t-statistics",col="lightblue")
hist(-t_twosample,main="t-statistics with pooled var",breaks=15,xlab="t-statistics with pooled variance",col="lightblue")

head(data.frame(t_stat=-t_twosample,moderated=final[,3]),n=10)

head(data.frame(p_value=p_twosample,moderated_pvalue=final[,4]),n=10)

Bonferroni-Holm-FDR

Bonferroni-Holm-FDR are applied first to original data and then to the transfromed data . This is done to have a better understanding of the selection of the most expressed genes and how the transformation of the data brings increase in the number of the expressed genes. Observe the Ven Diagram before and after.

rm(list = ls())
set.seed(100)
all_reject <- NULL
for (j in 1:50) {
  reject <- NULL
  for (i in 1:20000) {
    x <- rnorm(100, 0, 1)
    y <- rnorm(100, 0, 1)
    temp <- t.test(x, y, alternative = "two.sided")
    if (temp$p.value < 0.05)
      reject[i] <- 1
    else
      reject[i] <- 0
  }
  all_reject[j] <- sum(reject)
  print(j)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50

all_reject

 [1] 1016 1008 1022 1006  961 1058  953 1022 1016 1050 1032 1017 1023  990  972
[16] 1021 1023 1045 1005 1004 1060  997  993  974 1011  962  973 1017 1016 1024
[31] 1014 1044  937 1006 1038  995 1039 1002 1072 1034  974  993 1027 1007  997
[46] 1061 1031  995  966 1011

####Back to Colon example
library(plsgenomics)

For any news related to the 'plsgenomics' package (update, corrected bugs), please check http://thoth.inrialpes.fr/people/gdurif/
C++ based sparse PLS routines will soon be available on the CRAN in the new 'fastPLS' package.

data("Colon")

head(Colon$X[, 1:4], n = 10)

          1        2        3        4
1  8589.416 5468.241 4263.408 4064.936
2  9164.254 6719.529 4883.449 3718.159
3  3825.705 6970.361 5369.969 4705.650
4  6246.449 7823.534 5955.835 3975.564
5  3230.329 3694.450 3400.740 3463.586
6  2510.325 1960.655 1566.315 3072.816
7  7126.599 3779.068 3705.554 6594.514
8  4028.710 3156.159 2870.255 4417.591
9  9330.679 7017.230 4723.783 9491.534
10 5271.517 4740.768 3318.514 6792.348

colnames(Colon$X) <- paste("Gene", 1:ncol(Colon$X))

library(limma)
library(Biobase)
my_set <- ExpressionSet(assayData = t(Colon$X), phenodata = Colon$Y)

design <- model.matrix(~ as.factor(Colon$Y))
colnames(design) <- c("Intercept", "Tumor-Normal")

mod1 <- lmFit(my_set, design = design)
mod2 <- eBayes(mod1)
final <- topTable(mod2, coef=2, adjust="bonferroni",sort="p",n=Inf) 
head(final)

final[,4][1:6]*nrow(t(Colon$X))

[1] 5.276829e-05 4.658572e-04 6.927444e-04 7.498211e-04 9.878074e-04 1.558762e-03

final <- topTable(mod2, coef=2, adjust="holm",sort="p",n=Inf) 
head(final)

rank<-1:6
final[,4][1:6]*(nrow(t(Colon$X))-rank+1)

[1] 5.276829e-05 4.656243e-04 6.920517e-04 7.486964e-04 9.858317e-04 1.554865e-03

final <- topTable(mod2, coef=2, adjust="fdr",sort="p",n=Inf) 
head(final)

rank<-1:6
final[,4][1:6]*(nrow(t(Colon$X))/rank)

[1] 5.276829e-05 2.329286e-04 2.309148e-04 1.874553e-04 1.975615e-04 2.597936e-04

b1<-decideTests(mod2[,"Tumor-Normal"],adjust="bonferroni",p.value=0.05)
summary(b1)

       Tumor-Normal
Down              8
NotSig         1992
Up                0

b2<-decideTests(mod2[,"Tumor-Normal"],adjust="holm",p.value=0.05)
summary(b2)

       Tumor-Normal
Down              8
NotSig         1992
Up                0

b3<-decideTests(mod2[,"Tumor-Normal"],adjust="fdr",p.value=0.05)
summary(b3)

       Tumor-Normal
Down             24
NotSig         1937
Up               39

# So again here we observe that FDR is less conservative and because the other methods give just 8 genes we will decide to go with FDR decision of considering as differentially expressed genes 63.

ball<-b1
mat<-cbind(b1@.Data,b2@.Data,b3@.Data)
colnames(mat)<-c("Bonferroni","Holm","FDR")
ball@.Data<-mat

vennDiagram(b1,circle.col=c("orange"),lwd=2)

vennDiagram(b2,circle.col=c("orange"),lwd=2)

vennDiagram(b3,circle.col=c("orange"),lwd=2)

head(b3)

TestResults matrix
       Tumor-Normal
Gene 1            0
Gene 2            0
Gene 3            0
Gene 4            0
Gene 5            0
Gene 6            0

gene_ID<-which(b3==1)
gene_ID

 [1]   26   31   43   47   62   72   75  127  138  187  241  365  399  467  513
[16]  515  590  625  652  780  802  964  992 1002 1042 1047 1060 1067 1153 1293
[31] 1325 1346 1582 1634 1771 1772 1870 1972 1993

length(which(b3==1))# gives the position of genes which have value 1 corresponding to differentially expressed genes among Tumor and Normal (sample 1 and sample 2)

[1] 39

length(which(b3==0)) # # gives the position of genes which have value 0 corresponding to NOT differentially expressed gene

[1] 1937

par(mfrow=c(1,3))

vennDiagram(b1,circle.col=c("orange"),lwd=2, main="Bonferroni")
vennDiagram(b2,circle.col=c("orange"),lwd=2,main="Holm")

vennDiagram(b3,circle.col=c("orange"),lwd=2,main="FDR")

How we increased the number of differentially expressed genes by transforming the data? Observe above the set of genes when directly applied to the original data is less than the number of genes obtained after we use log2 transfom .

Log2 transformed data

data("Colon")
# Recall our notations:
expression.mat<-t(Colon$X)
log2_col.matrix <- log2(expression.mat)
# as.matrix(t(log2_col.matrix)) log2 matrix transformed

####Back to Colon example
library(plsgenomics)
data("Colon")

head(Colon$X[, 1:4], n = 10)

          1        2        3        4
1  8589.416 5468.241 4263.408 4064.936
2  9164.254 6719.529 4883.449 3718.159
3  3825.705 6970.361 5369.969 4705.650
4  6246.449 7823.534 5955.835 3975.564
5  3230.329 3694.450 3400.740 3463.586
6  2510.325 1960.655 1566.315 3072.816
7  7126.599 3779.068 3705.554 6594.514
8  4028.710 3156.159 2870.255 4417.591
9  9330.679 7017.230 4723.783 9491.534
10 5271.517 4740.768 3318.514 6792.348

colnames(Colon$X) <- paste("Gene", 1:ncol(Colon$X))

library(limma)
library(Biobase)
my_set <- ExpressionSet(assayData = log2_col.matrix, phenodata = Colon$Y)

design <- model.matrix(~ as.factor(Colon$Y))
colnames(design) <- c("Intercept", "Tumor-Normal")

mod1 <- lmFit(my_set, design = design)
mod2 <- eBayes(mod1)
final <- topTable(mod2, coef=2, adjust="bonferroni",sort="p",n=Inf) 
head(final)

final[,4][1:6]*nrow(t(Colon$X))

[1] 3.035236e-05 1.378366e-04 1.676019e-04 3.285367e-04 1.027929e-03 1.367052e-03

final <- topTable(mod2, coef=2, adjust="holm",sort="p",n=Inf) 
head(final)

rank<-1:6
final[,4][1:6]*(nrow(t(Colon$X))-rank+1)

[1] 3.035236e-05 1.377677e-04 1.674343e-04 3.280439e-04 1.025873e-03 1.363634e-03

final <- topTable(mod2, coef=2, adjust="fdr",sort="p",n=Inf) 
head(final)

rank<-1:6
final[,4][1:6]*(nrow(t(Colon$X))/rank)

[1] 3.035236e-05 6.891831e-05 5.586730e-05 8.213418e-05 2.055858e-04 2.278420e-04

b1<-decideTests(mod2[,"Tumor-Normal"],adjust="bonferroni",p.value=0.05)
summary(b1)

       Tumor-Normal
Down             12
NotSig         1976
Up               12

b2<-decideTests(mod2[,"Tumor-Normal"],adjust="holm",p.value=0.05)
summary(b2)

       Tumor-Normal
Down             12
NotSig         1975
Up               13

b3<-decideTests(mod2[,"Tumor-Normal"],adjust="fdr",p.value=0.05)
summary(b3)

       Tumor-Normal
Down             29
NotSig         1857
Up              114

# So again here we observe that FDR is less conservative and because the other methods give just 24,25 genes we will decide to go with FDR decision of considering as differentially expressed genes 143.

ball<-b1
mat<-cbind(b1@.Data,b2@.Data,b3@.Data)
colnames(mat)<-c("Bonferroni","Holm","FDR")
ball@.Data<-mat


par(mfrow=c(1,3))
vennDiagram(b1,circle.col=c("orange"),lwd=2, main="Bonferroni")
vennDiagram(b2,circle.col=c("orange"),lwd=2,main="Holm")

vennDiagram(b3,circle.col=c("orange"),lwd=2,main="FDR")

head(b3)

TestResults matrix
  Tumor-Normal
1            0
2            0
3            0
4            0
5            0
6            0

table(b3)

b3
  -1    0    1 
  29 1857  114 

gene_ID<-which(b3 %in% c(-1,1))
head(gene_ID)

[1] 26 31 43 47 62 66

length(gene_ID)# gives the position of genes which have value -1 or  1 corresponding to ( up /down) differentially expressed genes among Tumor and Normal (sample 1 and sample 2)

[1] 143

length(which(b3==0)) # # gives the position of genes which have value 0 corresponding to NOT differentially expressed gene

[1] 1857

Volcano Plot

NOTE IMPORTANT: Again here we have used the log2 transformed data for further processing. (not the original data)

Reference: https://combine-australia.github.io/RNAseq-R/06-rnaseq-day1.html#Plots_after_testing_for_DE

# We want to highlight the significant genes. We can get this from decideTests.
par(mfrow=c(1,3))
plotMD(mod2,coef=1,status=b1, values = c(-1, 1), hl.col=c("blue","red"),main="Bonferroni")
plotMD(mod2,coef=1,status=b2, values = c(-1, 1), hl.col=c("blue","red"),main="Holm")

plotMD(mod2,coef=1,status=b3, values = c(-1, 1), hl.col=c("blue","red"),main="FDR")

# For the volcano plot we have to specify how many of the top genes to highlight.
# We can also specify that we want to plot the gene symbol for the highlighted genes.
# let's highlight the top 100 most DE genes
par(mfrow=c(1,3))

volcanoplot(mod2,coef=1,highlight=100,names=mod2$genes$SYMBOL, main="....")
volcanoplot(mod2,coef=1,highlight=143,names=mod2$genes$SYMBOL, main="....")

volcanoplot(mod2,coef=1,highlight=300,names=mod2$genes$SYMBOL, main="....")

Question 2

Conduct PCA using differentially expressed genes identified from (1)

Now we will create a subset of the selected differentially expressed genes from FDR.

Colon_subset<-log2_col.matrix[gene_ID,]
View(Colon_subset)
dim(Colon_subset)# new dataset dimensions

[1] 143  62

Colon_sub_Y<-Colon$Y[gene_ID]
length(Colon_sub_Y) # the subset of genes from Colon$Y, length 114 genes

[1] 143

library(ComplexHeatmap)

Loading required package: grid
========================================
ComplexHeatmap version 2.12.1
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite either one:
- Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
    genomic data. Bioinformatics 2016.
- Gu, Z. Complex Heatmap Visualization. iMeta 2022.


The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================

Heatmap(as.matrix(Colon_subset),row_split = factor(Colon$Y[gene_ID]),cluster_columns = ,border=T,heatmap_legend_param = list(title = ""), cluster_row_slices = F)

library(plsgenomics)
library(Biobase)
library(limma)
data("Colon")
# Recall our notations:
expression.mat<-t(Colon$X)
log2_col.matrix <- log2(expression.mat)
# as.matrix(t(log2_col.matrix)) log2 matrix transformed

Reference: https://rpubs.com/namitakadam28/335503

To see which genes differ the most across the two groups(healthy,diseased), I perform the PCA on the data using the prcomp() function

By default, the prcomp() function centers the variables to have mean zero. By using the option scale=TRUE, we scale the variables to have standard deviation one.

The rotation matrix provides the principal component loadings; each column of pr.out$rotation contains the corresponding principal component loading vector.

In this case, the loading can be considered to be the weight of each gene in both the groups.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.8     ✔ stringr 1.4.1
✔ tidyr   1.2.1     ✔ forcats 0.5.2
✔ readr   2.1.3     ── Conflicts ──────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::combine()    masks Biobase::combine(), BiocGenerics::combine()
✖ dplyr::filter()     masks stats::filter()
✖ dplyr::lag()        masks stats::lag()
✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()

pr.out <- prcomp(log2_col.matrix[gene_ID,], center=TRUE,scale = TRUE)
head(pr.out$rotation)

         PC1         PC2        PC3        PC4          PC5         PC6
1 -0.1371324  0.01424307 -0.1881907 0.14980265 -0.093092513  0.03414976
2 -0.1197226 -0.17269200 -0.2081417 0.09746436  0.006173545 -0.03202835
3 -0.1181497 -0.06364682 -0.1030101 0.08218244 -0.222324083  0.43884363
4 -0.1242051 -0.10582203 -0.1758494 0.13677665 -0.177866248  0.10260358
5 -0.1376373  0.04513975  0.1384985 0.01233534 -0.080547646 -0.07462775
6 -0.1188489 -0.16381324  0.2140853 0.12292452 -0.065666815 -0.04828333
          PC7         PC8          PC9         PC10        PC11        PC12
1 -0.13701632 -0.03515201  0.088875021 -0.053618097  0.03744552 -0.21954688
2 -0.05352914 -0.01377077  0.048234206 -0.033501555  0.06322502 -0.16824845
3  0.03049032  0.06270598 -0.174917800  0.098012155 -0.32299011  0.03149722
4 -0.24073543  0.06833523 -0.082352135  0.076937232 -0.05554723 -0.01023646
5  0.24039817  0.07255428  0.005720992 -0.104506373  0.02632018 -0.19873993
6 -0.02427577  0.04329205  0.110166766  0.007395204  0.03285915 -0.05181883
         PC13        PC14        PC15        PC16        PC17         PC18
1  0.11880784 -0.03192956  0.06169239 -0.16672404  0.06448241 -0.110203300
2 -0.04374693  0.13493530 -0.15993067  0.13827454 -0.06460728 -0.005700702
3  0.07145012 -0.25587057  0.07433530 -0.06744187  0.10705375 -0.138844801
4 -0.18571840  0.02748144 -0.05601846  0.05497411 -0.13936323  0.022923233
5 -0.07186679 -0.28032161 -0.19012318  0.04318114 -0.11272248  0.140389635
6  0.12966516 -0.07113824 -0.03248420  0.08041356  0.08504552  0.101941352
         PC19          PC20        PC21        PC22        PC23        PC24
1  0.11607997 -0.0756183391 -0.07843758  0.23330602  0.04121391  0.01183369
2 -0.09355258 -0.0500788818 -0.03750728  0.01392347 -0.05235914 -0.15706909
3  0.01612754 -0.0651637760  0.23871092 -0.04493530 -0.01424829 -0.16770155
4 -0.18413531  0.1013959597 -0.00783001 -0.45426157  0.16353780  0.03707161
5 -0.07520814  0.2053277177 -0.08398741 -0.12329769  0.17394965 -0.05617818
6  0.10473191  0.0002324725 -0.02040530 -0.01518136 -0.14522440 -0.14265632
         PC25        PC26        PC27         PC28        PC29        PC30
1 -0.22532769 -0.09620743  0.14481715  0.061328389  0.22456470 -0.01354592
2  0.07912227 -0.01866530  0.08501488  0.005927554  0.15720180 -0.06473936
3  0.19632050 -0.02821033  0.02605289  0.068797302  0.12720073  0.30990693
4  0.07299823  0.08449219  0.21435717 -0.102401173 -0.21126478  0.02976040
5 -0.06791550  0.04571789  0.05460287  0.188647985  0.23475370 -0.01727302
6 -0.14805882  0.08560316 -0.13705693  0.061280949  0.07093735 -0.10152201
         PC31         PC32        PC33        PC34        PC35         PC36
1 -0.12988685 -0.006223575 -0.25660025  0.03287804 0.126896725  0.216797180
2 -0.09261889  0.131428461  0.02141692  0.12074323 0.009273666 -0.022834927
3  0.04198371  0.099610465  0.09780292 -0.11950939 0.115412140 -0.063202993
4  0.05902507  0.057424206  0.07582962  0.00379116 0.033286002 -0.149765300
5 -0.13049547 -0.179884166  0.01741862  0.10271220 0.050578712 -0.008278755
6 -0.14871141  0.017599205  0.08282409 -0.51452499 0.036453373 -0.153791755
         PC37         PC38        PC39       PC40        PC41        PC42
1  0.07577961 -0.179353281 -0.03512746 -0.2026418  0.10133392  0.08862831
2 -0.15085024 -0.090259648  0.14483226 -0.1245761  0.08635337  0.07203724
3 -0.05314329  0.085799991  0.10473527  0.1283999 -0.05309397  0.09840348
4  0.27158705 -0.086096579  0.05159156 -0.2019021  0.07476713 -0.07611619
5  0.39786870 -0.008727689 -0.04170150  0.2343805 -0.02745676  0.10670243
6  0.04923080  0.252949085  0.08310211 -0.2112687  0.12391785  0.02821809
         PC43         PC44        PC45        PC46        PC47        PC48
1  0.05767448  0.195335733 -0.01199872 -0.06701675 -0.05229611 -0.05831927
2  0.09340344  0.022121642 -0.07172321  0.14743996 -0.03840555 -0.08193686
3  0.10092785  0.137151195 -0.06118984  0.08178975 -0.12269027 -0.09631932
4 -0.11530163  0.004110503  0.04495368 -0.09458246  0.09777637  0.02612710
5  0.20681916  0.045137687  0.07616930  0.06284382 -0.01222206  0.10890470
6 -0.12955897 -0.033042690  0.23877925  0.03663674 -0.13406934  0.10003305
         PC49        PC50        PC51         PC52         PC53         PC54
1 -0.09181136 -0.03228224  0.03512305  0.007548056  0.301907573 -0.184352760
2  0.02259207  0.13490503  0.14529173  0.125351413  0.001973985  0.313140313
3  0.04222169 -0.04472071  0.03946985  0.128622932 -0.047344866 -0.068606707
4  0.01468639  0.07230853 -0.05873750 -0.168312783  0.090730674  0.003430529
5 -0.06736195 -0.04578512 -0.17968898  0.156007087 -0.051815576  0.130582690
6  0.08608957  0.02393487  0.06186330 -0.064117976 -0.068958067 -0.106045599
         PC55         PC56         PC57         PC58        PC59         PC60
1 -0.04077794 -0.130921377 -0.064899235 -0.173555029  0.17041967  0.048501507
2  0.10199776 -0.151944589  0.150677223  0.386706239  0.01450419 -0.355484884
3 -0.05574786 -0.021053598 -0.029855070 -0.008982396 -0.03543677  0.007789275
4  0.09327982 -0.123903465  0.042824902 -0.075291602  0.05813638  0.078634420
5  0.03032557  0.130998949 -0.017590328 -0.033233264  0.09994172 -0.066835064
6  0.21889393 -0.009786299 -0.008418231 -0.108699319 -0.02048189 -0.037518150
         PC61        PC62
1  0.18530379 -0.02446595
2 -0.08434292  0.22853130
3 -0.10076178 -0.02149227
4  0.22772838 -0.05946191
5 -0.03726423 -0.01121972
6  0.07512160  0.22185984

# apply() function allows us to apply a function-in this case, the sum() function-to each row of the data set
# The second input here denotes whether we wish to compute the sum of the rows, 1, or the columns, 2.
# So, I am summing all the weights of each gene to know what is the total weightage of each gene in this distribution

gene.load <- apply(pr.out$rotation, 1, sum)

# Now, I am arranging the weightage obtained above in a descending order by taking the absolute value of the total loadings for each gene
gene.differ <- order(abs(gene.load), decreasing = TRUE)

# To show the top most different genes across the two groups
gene.differ[1:50]

 [1] 10 34 51 13 16 24 25  7 31 62 32 11 61 19 23 43 60 26 48  5 59 33 47 30 20  4
[27] 17 18 15 56 54 57 55  9  3  2 45 42 22 28 27 52 40 21  1 44 53 41 46 37

#display principal components
biplot(pr.out)

screeplot(pr.out)

Reference: http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/

library("FactoMineR")

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

res.pca <- PCA(log2_col.matrix[gene_ID,], graph = FALSE)
print(res.pca)

**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 143 individuals, described by 62 variables
*The results are available in the following objects:

   name               description                          
1  "$eig"             "eigenvalues"                        
2  "$var"             "results for the variables"          
3  "$var$coord"       "coord. for the variables"           
4  "$var$cor"         "correlations variables - dimensions"
5  "$var$cos2"        "cos2 for the variables"             
6  "$var$contrib"     "contributions of the variables"     
7  "$ind"             "results for the individuals"        
8  "$ind$coord"       "coord. for the individuals"         
9  "$ind$cos2"        "cos2 for the individuals"           
10 "$ind$contrib"     "contributions of the individuals"   
11 "$call"            "summary statistics"                 
12 "$call$centre"     "mean of the variables"              
13 "$call$ecart.type" "standard error of the variables"    
14 "$call$row.w"      "weights for the individuals"        
15 "$call$col.w"      "weights for the variables"          

library("factoextra")

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

eig.val <- get_eigenvalue(res.pca)
eig.val

         eigenvalue variance.percent cumulative.variance.percent
Dim.1  41.610771938     67.114148286                    67.11415
Dim.2   7.983353186     12.876376107                    79.99052
Dim.3   2.197108892      3.543724019                    83.53425
Dim.4   1.485589704      2.396112425                    85.93036
Dim.5   1.250328032      2.016658117                    87.94702
Dim.6   0.907472116      1.463664703                    89.41068
Dim.7   0.625516684      1.008897878                    90.41958
Dim.8   0.533680435      0.860774895                    91.28036
Dim.9   0.426978388      0.688674819                    91.96903
Dim.10  0.372913567      0.601473495                    92.57050
Dim.11  0.319103653      0.514683312                    93.08519
Dim.12  0.307742901      0.496359518                    93.58155
Dim.13  0.283892622      0.457891326                    94.03944
Dim.14  0.273865270      0.441718177                    94.48116
Dim.15  0.242092197      0.390471286                    94.87163
Dim.16  0.213792953      0.344827344                    95.21646
Dim.17  0.192121125      0.309872782                    95.52633
Dim.18  0.176796134      0.285155055                    95.81148
Dim.19  0.165401563      0.266776715                    96.07826
Dim.20  0.152381354      0.245776377                    96.32404
Dim.21  0.149276729      0.240768918                    96.56481
Dim.22  0.140054280      0.225894001                    96.79070
Dim.23  0.132561842      0.213809422                    97.00451
Dim.24  0.125331257      0.202147189                    97.20666
Dim.25  0.114294731      0.184346341                    97.39100
Dim.26  0.107967025      0.174140363                    97.56514
Dim.27  0.100995734      0.162896346                    97.72804
Dim.28  0.096940096      0.156354994                    97.88439
Dim.29  0.091413044      0.147440394                    98.03183
Dim.30  0.084984740      0.137072162                    98.16891
Dim.31  0.080365499      0.129621773                    98.29853
Dim.32  0.077626500      0.125204032                    98.42373
Dim.33  0.076489025      0.123369395                    98.54710
Dim.34  0.069672497      0.112374996                    98.65948
Dim.35  0.066232197      0.106826124                    98.76630
Dim.36  0.061205218      0.098718094                    98.86502
Dim.37  0.060133609      0.096989692                    98.96201
Dim.38  0.055316207      0.089219688                    99.05123
Dim.39  0.053210292      0.085823051                    99.13705
Dim.40  0.045882055      0.074003315                    99.21106
Dim.41  0.043519175      0.070192218                    99.28125
Dim.42  0.040590299      0.065468224                    99.34672
Dim.43  0.039751744      0.064115716                    99.41083
Dim.44  0.035421802      0.057131939                    99.46797
Dim.45  0.034368484      0.055433039                    99.52340
Dim.46  0.030309114      0.048885668                    99.57228
Dim.47  0.027299088      0.044030787                    99.61631
Dim.48  0.026221117      0.042292124                    99.65861
Dim.49  0.025008309      0.040335983                    99.69894
Dim.50  0.023929793      0.038596440                    99.73754
Dim.51  0.020831718      0.033599545                    99.77114
Dim.52  0.019458895      0.031385315                    99.80252
Dim.53  0.018681166      0.030130913                    99.83265
Dim.54  0.016342746      0.026359267                    99.85901
Dim.55  0.015494985      0.024991911                    99.88401
Dim.56  0.014028664      0.022626878                    99.90663
Dim.57  0.012218660      0.019707516                    99.92634
Dim.58  0.011756405      0.018961944                    99.94530
Dim.59  0.010528501      0.016981453                    99.96228
Dim.60  0.010036363      0.016187683                    99.97847
Dim.61  0.007162012      0.011551633                    99.99002
Dim.62  0.006185665      0.009976879                   100.00000

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 85))

var <- get_pca_var(res.pca)
var

Principal Component Analysis Results for variables
 ===================================================
  Name       Description                                    
1 "$coord"   "Coordinates for the variables"                
2 "$cor"     "Correlations between variables and dimensions"
3 "$cos2"    "Cos2 for the variables"                       
4 "$contrib" "contributions of the variables"               

The correlation between a variable and a principal component (PC) is used as the coordinates of the variable on the PC. The representation of variables differs from the plot of the observations: The observations are represented by their projections, but the variables are represented by their correlations (Abdi and Williams 2010).

# Coordinates of variables
head(var$coord, 4)

      Dim.1       Dim.2      Dim.3      Dim.4        Dim.5
1 0.8845917 -0.04024355 -0.2789484 -0.1825866 -0.104094249
2 0.7722874  0.48793828 -0.3085211 -0.1187942  0.006903139
3 0.7621416  0.17983298 -0.1526883 -0.1001679 -0.248598494
4 0.8012028  0.29899834 -0.2606553 -0.1667099 -0.198886602

Quality of representation

The quality of representation of the variables on factor map is called cos2 (square cosine, squared coordinates) . You can access to the cos2 as follow:

head(var$cos2, 4)

      Dim.1       Dim.2      Dim.3      Dim.4        Dim.5
1 0.7825025 0.001619544 0.07781223 0.03333787 1.083561e-02
2 0.5964278 0.238083766 0.09518525 0.01411206 4.765333e-05
3 0.5808598 0.032339902 0.02331370 0.01003360 6.180121e-02
4 0.6419259 0.089400006 0.06794121 0.02779219 3.955588e-02

# Total cos2 of variables on Dim.1 and Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)

The plot below is also known as variable correlation plots. It shows the relationships between all variables. It can be interpreted as follow:

Positively correlated variables are grouped together. Negatively correlated variables are positioned on opposite sides of the plot origin (opposed quadrants). The distance between variables and the origin measures the quality of the variables on the factor map. Variables that are away from the origin are well represented on the factor map.

Note that,

A high cos2 indicates a good representation of the variable on the principal component. In this case the variable is positioned close to the circumference of the correlation circle.

A low cos2 indicates that the variable is not perfectly represented by the PCs. In this case the variable is close to the center of the circle.

It’s possible to color variables by their cos2 values using the argument col.var = “cos2”. This produces a gradient colors. In this case, the argument gradient.cols can be used to provide a custom color. For instance, gradient.cols = c(“white”, “blue”, “red”) means that:

variables with low cos2 values will be colored in “white” variables with mid cos2 values will be colored in “blue” variables with high cos2 values will be colored in red

# Color by cos2 values: quality on the factor map
fviz_pca_var(res.pca, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE # Avoid text overlapping
             )

Contributions of variables to PCs

The contributions of variables in accounting for the variability in a given principal component are expressed in percentage.

Variables that are correlated with PC1 (i.e., Dim.1) and PC2 (i.e., Dim.2) are the most important in explaining the variability in the data set. Variables that do not correlated with any PC or correlated with the last dimensions are variables with low contribution and might be removed to simplify the overall analysis. The larger the value of the contribution, the more the variable contributes to the component. The contribution of variables can be extracted as follow :

head(var$contrib, 4)

     Dim.1      Dim.2    Dim.3     Dim.4       Dim.5
1 1.880529 0.02028651 3.541573 2.2440835 0.866621590
2 1.433350 2.98225270 4.332296 0.9499302 0.003811266
3 1.395936 0.40509172 1.061108 0.6753954 4.942799788
4 1.542692 1.11983028 3.092301 1.8707852 3.163640229

library("corrplot")
corrplot(var$contrib[1:5,], is.corr=FALSE) 


par(mfrow=c(1,6))

 corrplot(var$contrib[1:10,], is.corr=FALSE) 
 corrplot(var$contrib[11:20,], is.corr=FALSE) 

 corrplot(var$contrib[21:30,], is.corr=FALSE) 
 corrplot(var$contrib[31:40,], is.corr=FALSE) 

 corrplot(var$contrib[41:50,], is.corr=FALSE) 
 corrplot(var$contrib[51:62,], is.corr=FALSE)

 
 
 
 par(mfrow=c(1,6))

 corrplot(var$cor[1:10,], is.corr=FALSE) 
 corrplot(var$cor[11:20,], is.corr=FALSE) 

 corrplot(var$cor[21:30,], is.corr=FALSE) 
 corrplot(var$cor[31:40,], is.corr=FALSE) 

 corrplot(var$cor[41:50,], is.corr=FALSE) 
 corrplot(var$cor[51:62,], is.corr=FALSE)

Below the contribution of variables to dim 1 and dim 2 and both dim are shown by a bar chart.

# Contributions of variables to PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)

fviz_contrib(res.pca, choice = "var", axes = 1)

# Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)

# Total contribution PC1 and PC2
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)

NA
NA

The most important (or, contributing) variables can be highlighted on the correlation plot as follow:

fviz_pca_var(res.pca, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
             )

# Create a random continuous variable of length 10
set.seed(123)
my.cont.var <- rnorm(62)
# Color variables by the continuous variable
fviz_pca_var(res.pca, col.var = my.cont.var,
             gradient.cols = c("blue", "yellow", "red"),
             legend.title = "Cont.Var")

# Create a grouping variable using kmeans
# Create 2 groups of variables (centers = 2)
set.seed(123)
res.km <- kmeans(var$coord, centers = 2, nstart = 25)
grp <- as.factor(res.km$cluster)
# Color variables by groups
fviz_pca_var(res.pca, col.var = grp,  palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),legend.title = "Cluster")

# Create a grouping variable using kmeans
# Create 3 groups of variables (centers = 3)
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)
# Color variables by groups
fviz_pca_var(res.pca, col.var = grp,   palette = c("#0073C2FF", "#EFC000FF", "#868686FF"), legend.title = "Cluster")

Dimension description

Above we described how to highlight variables according to their contributions to the principal components.

Note also that, the function dimdesc() [in FactoMineR], for dimension description, can be used to identify the most significantly associated variables with a given principal component . It can be used as follow:

res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
# Description of dimension 1
res.desc$Dim.1

Link between the variable and the continuous variables (R-square)
=================================================================================
   correlation      p.value
13   0.9263233 1.200862e-61
15   0.9122200 1.673564e-56
34   0.9085283 2.676463e-55
53   0.9029555 1.417047e-53
35   0.8880039 2.015341e-49
5    0.8878490 2.209094e-49
31   0.8857127 7.732266e-49
27   0.8848445 1.277469e-48
1    0.8845917 1.477437e-48
25   0.8727034 9.624798e-46
41   0.8719548 1.415931e-45
32   0.8669709 1.741780e-44
33   0.8652106 4.125369e-44
52   0.8627879 1.325280e-43
51   0.8594327 6.433540e-43
40   0.8590668 7.624106e-43
17   0.8547057 5.566024e-42
7    0.8518120 2.008971e-41
14   0.8476333 1.222736e-40
45   0.8473754 1.364475e-40
59   0.8461740 2.268388e-40
16   0.8435396 6.811561e-40
30   0.8401878 2.680440e-39
46   0.8344248 2.628829e-38
58   0.8332144 4.199406e-38
24   0.8323952 5.753629e-38
18   0.8317913 7.248864e-38
28   0.8272517 3.999427e-37
44   0.8269230 4.517037e-37
61   0.8259456 6.477181e-37
21   0.8145028 3.752594e-35
12   0.8139546 4.526111e-35
62   0.8121873 8.246768e-35
29   0.8075920 3.811283e-34
4    0.8012028 2.993860e-33
37   0.8000633 4.290028e-33
38   0.7977045 8.967425e-33
8    0.7975281 9.472178e-33
19   0.7954779 1.782670e-32
50   0.7903535 8.394472e-32
47   0.7892055 1.180797e-31
36   0.7871540 2.161176e-31
42   0.7831828 6.834730e-31
43   0.7819325 9.771930e-31
10   0.7803765 1.519735e-30
26   0.7790859 2.186093e-30
20   0.7761213 4.993053e-30
2    0.7722874 1.426246e-29
22   0.7710343 2.001016e-29
56   0.7704433 2.345785e-29
55   0.7687748 3.665066e-29
6    0.7666514 6.432215e-29
3    0.7621416 2.083017e-28
54   0.7578049 6.293198e-28
9    0.7564266 8.899658e-28
23   0.7559212 1.009969e-27
57   0.7495143 4.890215e-27
49   0.7449073 1.476452e-26
39   0.7441826 1.752966e-26
60   0.7420981 2.862995e-26
11   0.7337655 1.942641e-25
48   0.7080511 4.666297e-23

#Graph of individuals

ind <- get_pca_ind(res.pca)
ind

Principal Component Analysis Results for individuals
 ===================================================
  Name       Description                       
1 "$coord"   "Coordinates for the individuals" 
2 "$cos2"    "Cos2 for the individuals"        
3 "$contrib" "contributions of the individuals"

# Coordinates of individuals
head(ind$coord)

      Dim.1      Dim.2     Dim.3       Dim.4      Dim.5
26 18.90011 -1.0919715 0.3039404  0.41098891 -0.3376709
31 15.28620 -1.2149321 1.1250378 -0.24045752 -0.2921266
43 13.99193 -0.9685592 1.2170596 -0.17396722 -0.4465271
47 14.85878 -1.1133445 1.1399760 -0.78381212 -0.9176131
62 12.68649 -1.1948691 1.1984574  0.11631039 -0.3640826
66 10.44386  3.7588282 1.5735257  0.05247403 -0.3290574

# Quality of individuals
head(ind$cos2)

       Dim.1       Dim.2        Dim.3        Dim.4        Dim.5
26 0.9911947 0.003308665 0.0002563346 4.686951e-04 0.0003163862
31 0.9809147 0.006196357 0.0053133280 2.427220e-04 0.0003582406
43 0.9809283 0.004700397 0.0074217423 1.516411e-04 0.0009990284
47 0.9687967 0.005439081 0.0057024015 2.695817e-03 0.0036947544
62 0.9707691 0.008611391 0.0086631893 8.159615e-05 0.0007995258
66 0.8343652 0.108078369 0.0189400669 2.106312e-05 0.0008282800

# Contributions of individuals
head(ind$contrib)

      Dim.1      Dim.2      Dim.3       Dim.4      Dim.5
26 6.003258 0.10444826 0.02940285 0.079510646 0.06377161
31 3.926965 0.12929527 0.40285299 0.027217098 0.04772900
43 3.290131 0.08217335 0.47145045 0.014246229 0.11151560
47 3.710432 0.10857698 0.41362215 0.289193833 0.47093306
62 2.704839 0.12506026 0.45714876 0.006367986 0.07413783
66 1.833076 1.23760933 0.78806141 0.001296146 0.06055965

Plots: quality and contribution

The fviz_pca_ind() is used to produce the graph of individuals. To create a simple plot, type this:

fviz_pca_ind(res.pca)

#Like variables, it’s also possible to color individuals by their cos2 values:

fviz_pca_ind(res.pca, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

NOTE: that, individuals that are similar are grouped together on the plot.

fviz_pca_ind(res.pca, pointsize = "cos2", 
             pointshape = 21, fill = "#E7B800",
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

To create a bar plot of the quality of representation (cos2) of individuals on the factor map, you can use the function fviz_cos2() as previously described for variables:

fviz_cos2(res.pca, choice = "ind")

To visualize the contribution of individuals to the first two principal components, type this:

# Total contribution on PC1 and PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", # Variables color
                col.ind = "#696969"  # Individuals color
                )

Question 3

Cluster Hierarchial

dist<-dist(log2_col.matrix[gene_ID,])

Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘spam’

trial_hclust<-hclust(dist)
plot(trial_hclust)

plot(trial_hclust,labels=Colon$gene.names[gene_ID])
abline(h = 20, col = "brown", lwd = 2) # add horizontal line to illustrate cutting dendrogram

pred_class<-cutree(trial_hclust,2)
pred_class

  26   31   43   47   62   66   72   75  100  102  127  138  141  143  147  187 
   1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
 190  201  241  245  249  264  267  281  286  295  343  365  377  391  399  415 
   1    1    1    1    1    1    1    1    1    1    1    1    1    1    2    1 
 444  489  493  495  513  515  520  529  549  550  571  581  590  601  619  625 
   1    1    1    1    1    2    1    1    1    2    2    1    1    1    2    1 
 639  648  652  698  739  765  780  802  822  830  882  897  964  989  992  994 
   2    1    1    2    1    1    1    1    1    1    1    1    2    2    1    2 
1002 1042 1047 1060 1067 1073 1110 1115 1153 1168 1187 1196 1208 1221 1227 1248 
   1    2    1    2    2    1    2    2    1    2    2    1    2    2    2    2 
1256 1258 1263 1293 1325 1334 1340 1346 1387 1406 1411 1414 1423 1466 1472 1473 
   2    1    2    1    1    2    2    2    1    2    1    1    2    2    2    2 
1489 1494 1511 1526 1546 1549 1582 1583 1634 1635 1637 1648 1659 1668 1671 1674 
   2    2    2    1    2    2    2    2    2    2    2    2    2    2    2    2 
1675 1679 1724 1730 1758 1763 1770 1771 1772 1785 1808 1836 1839 1843 1867 1870 
   2    2    2    2    2    2    2    2    2    2    2    1    2    2    2    2 
1884 1886 1892 1897 1900 1904 1920 1935 1959 1960 1967 1972 1974 1983 1993 
   2    2    2    2    2    2    2    2    2    2    2    2    2    2    2 

table(pred_class,Colon$gene.names[gene_ID])

          
pred_class Hsa.103 Hsa.1045 Hsa.1047 Hsa.10664 Hsa.1073 Hsa.1130 Hsa.1131
         1       1        1        1         0        1        1        1
         2       0        0        0         1        0        0        0
          
pred_class Hsa.1132 Hsa.1143 Hsa.1145 Hsa.11616 Hsa.11673 Hsa.1198 Hsa.1205
         1        1        0        0         0         1        0        1
         2        0        1        1         1         0        1        0
          
pred_class Hsa.1207 Hsa.1221 Hsa.1387 Hsa.14069 Hsa.1410 Hsa.1423 Hsa.1435
         1        1        1        1         0        0        0        0
         2        0        0        0         1        1        1        1
          
pred_class Hsa.1454 Hsa.1588 Hsa.1591 Hsa.1617 Hsa.1660 Hsa.168 Hsa.17564
         1        0        1        1        0        0       1         0
         2        1        0        0        1        1       0         1
          
pred_class Hsa.1832 Hsa.1902 Hsa.1985 Hsa.2097 Hsa.2250 Hsa.22762 Hsa.229
         1        0        1        1        0        0         0       0
         2        1        0        0        1        1         1       1
          
pred_class Hsa.2291 Hsa.2316 Hsa.2344 Hsa.2451 Hsa.2456 Hsa.24582 Hsa.25322
         1        0        0        0        0        0         0         0
         2        1        1        1        1        1         1         1
          
pred_class Hsa.2644 Hsa.2645 Hsa.26528 Hsa.2705 Hsa.2715 Hsa.27686 Hsa.2773
         1        0        0         1        0        1         1        1
         2        1        1         0        1        0         0        0
          
pred_class Hsa.2795 Hsa.2800 Hsa.2818 Hsa.2821 Hsa.2842 Hsa.286 Hsa.2863 Hsa.2928
         1        1        1        0        0        0       1        1        0
         2        0        0        1        1        1       0        0        1
          
pred_class Hsa.2950 Hsa.3001 Hsa.3007 Hsa.3016 Hsa.3083 Hsa.3088 Hsa.3093
         1        1        1        1        1        0        1        0
         2        0        0        0        0        1        0        1
          
pred_class Hsa.3152 Hsa.31630 Hsa.3252 Hsa.3263 Hsa.3296 Hsa.33 Hsa.3305 Hsa.3306
         1        1         1        0        0        0      0        1        1
         2        0         0        1        1        1      1        0        0
          
pred_class Hsa.3331 Hsa.3566 Hsa.36689 Hsa.36696 Hsa.36952 Hsa.37937 Hsa.3803
         1        0        1         1         0         1         1        1
         2        1        0         0         1         0         0        0
          
pred_class Hsa.3963 Hsa.39753 Hsa.40063 Hsa.41260 Hsa.41283 Hsa.41323 Hsa.4252
         1        0         0         1         0         0         0        1
         2        1         1         0         1         1         1        0
          
pred_class Hsa.43279 Hsa.451 Hsa.45446 Hsa.462 Hsa.466 Hsa.4689 Hsa.490 Hsa.491
         1         1       1         0       0       0        1       0       0
         2         0       0         1       1       1        0       1       1
          
pred_class Hsa.5211 Hsa.5392 Hsa.5398 Hsa.5444 Hsa.549 Hsa.558 Hsa.56 Hsa.5971
         1        1        0        1        1       0       1      0        0
         2        0        1        0        0       1       0      1        1
          
pred_class Hsa.6 Hsa.601 Hsa.6030 Hsa.627 Hsa.6458 Hsa.6472 Hsa.662 Hsa.678
         1     0       0        0       0        0        0       0       1
         2     1       1        1       1        1        1       1       0
          
pred_class Hsa.6782 Hsa.6814 Hsa.692 Hsa.7 Hsa.702 Hsa.726 Hsa.7395 Hsa.773
         1        0        0       3     0       1       0        1       1
         2        1        1       0     1       0       1        0       0
          
pred_class Hsa.8010 Hsa.8040 Hsa.8068 Hsa.810 Hsa.8125 Hsa.8147 Hsa.816 Hsa.821
         1        1        0        1       1        1        1       0       1
         2        0        1        0       0        0        0       1       0
          
pred_class Hsa.8214 Hsa.831 Hsa.832 Hsa.853 Hsa.878 Hsa.8781 Hsa.9218 Hsa.9353
         1        0       1       1       0       1        0        0        0
         2        1       0       0       1       0        1        1        1
          
pred_class Hsa.94 Hsa.951 Hsa.957 Hsa.960 Hsa.9744 Hsa.996 Hsa.9972
         1      0       1       1       0        0       1        0
         2      1       0       0       1        1       0        1

distance <- get_dist(log2_col.matrix[gene_ID,])

Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘spam’

fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high =
"#FC4E07"))

Reference: https://tavareshugo.github.io/data-carpentry-rnaseq/04b_rnaseq_clustering.html

gene_cluster <- cutree(trial_hclust, k = 5) %>% 
  # turn the named vector into a tibble
  enframe() 

head(gene_cluster)

table(gene_cluster)# will show in which cluster the gene is 

      value
name   1 2 3 4 5
  100  1 0 0 0 0
  1002 0 1 0 0 0
  102  1 0 0 0 0
  1042 0 0 1 0 0
  1047 0 1 0 0 0
  1060 0 0 1 0 0
  1067 0 0 1 0 0
  1073 0 1 0 0 0
  1110 0 0 1 0 0
  1115 0 0 1 0 0
  1153 0 1 0 0 0
  1168 0 0 1 0 0
  1187 0 0 1 0 0
  1196 0 1 0 0 0
  1208 0 0 0 1 0
  1221 0 0 1 0 0
  1227 0 0 1 0 0
  1248 0 0 1 0 0
  1256 0 0 1 0 0
  1258 0 1 0 0 0
  1263 0 0 1 0 0
  127  0 1 0 0 0
  1293 0 1 0 0 0
  1325 0 1 0 0 0
  1334 0 0 1 0 0
  1340 0 0 1 0 0
  1346 0 0 0 1 0
  138  1 0 0 0 0
  1387 0 1 0 0 0
  1406 0 0 1 0 0
  141  0 1 0 0 0
  1411 0 1 0 0 0
  1414 0 1 0 0 0
  1423 0 0 0 0 1
  143  0 1 0 0 0
  1466 0 0 0 1 0
  147  0 1 0 0 0
  1472 0 0 1 0 0
  1473 0 0 1 0 0
  1489 0 0 1 0 0
  1494 0 0 0 0 1
  1511 0 0 1 0 0
  1526 0 1 0 0 0
  1546 0 0 0 1 0
  1549 0 0 1 0 0
  1582 0 0 0 1 0
  1583 0 0 0 1 0
  1634 0 0 1 0 0
  1635 0 0 0 1 0
  1637 0 0 1 0 0
  1648 0 0 0 1 0
  1659 0 0 1 0 0
  1668 0 0 0 1 0
  1671 0 0 1 0 0
  1674 0 0 0 1 0
  1675 0 0 0 1 0
  1679 0 0 0 1 0
  1724 0 0 0 1 0
  1730 0 0 1 0 0
  1758 0 0 0 1 0
  1763 0 0 0 1 0
  1770 0 0 0 1 0
  1771 0 0 0 1 0
  1772 0 0 0 1 0
  1785 0 0 0 1 0
  1808 0 0 1 0 0
  1836 0 1 0 0 0
  1839 0 0 0 1 0
  1843 0 0 0 0 1
  1867 0 0 0 1 0
  187  1 0 0 0 0
  1870 0 0 1 0 0
  1884 0 0 0 1 0
  1886 0 0 0 1 0
  1892 0 0 0 1 0
  1897 0 0 0 1 0
  190  0 1 0 0 0
  1900 0 0 1 0 0
  1904 0 0 0 1 0
  1920 0 0 0 1 0
  1935 0 0 0 1 0
  1959 0 0 1 0 0
  1960 0 0 0 1 0
  1967 0 0 0 0 1
  1972 0 0 1 0 0
  1974 0 0 0 0 1
  1983 0 0 0 1 0
  1993 0 0 1 0 0
  201  0 1 0 0 0
  241  0 1 0 0 0
  245  0 1 0 0 0
  249  0 1 0 0 0
  26   1 0 0 0 0
  264  0 1 0 0 0
  267  0 1 0 0 0
  281  0 1 0 0 0
  286  1 0 0 0 0
  295  0 1 0 0 0
  31   1 0 0 0 0
  343  0 1 0 0 0
  365  0 1 0 0 0
  377  0 1 0 0 0
  391  0 1 0 0 0
  399  0 0 1 0 0
  415  0 1 0 0 0
  43   1 0 0 0 0
  444  0 1 0 0 0
  47   1 0 0 0 0
  489  0 1 0 0 0
  493  0 1 0 0 0
  495  0 1 0 0 0
  513  0 1 0 0 0
  515  0 0 1 0 0
  520  0 1 0 0 0
  529  0 1 0 0 0
  549  0 1 0 0 0
  550  0 0 1 0 0
  571  0 0 1 0 0
  581  0 1 0 0 0
  590  0 1 0 0 0
  601  0 1 0 0 0
  619  0 0 1 0 0
  62   1 0 0 0 0
  625  0 1 0 0 0
  639  0 0 1 0 0
  648  0 1 0 0 0
  652  0 1 0 0 0
  66   1 0 0 0 0
  698  0 0 1 0 0
  72   1 0 0 0 0
  739  0 1 0 0 0
  75   1 0 0 0 0
  765  0 1 0 0 0
  780  0 1 0 0 0
  802  0 1 0 0 0
  822  0 1 0 0 0
  830  0 1 0 0 0
  882  0 1 0 0 0
  897  0 1 0 0 0
  964  0 0 1 0 0
  989  0 0 1 0 0
  992  0 1 0 0 0
  994  0 0 1 0 0

plot(trial_hclust, cex = 0.6)
rect.hclust(trial_hclust, k = 4, border = 2:4)

plot(trial_hclust, cex = 0.6)
rect.hclust(trial_hclust, k = 3, border = 2:3)

#Visualise gene expression trends per cluster

trans_cts_cluster <- cbind(log2_col.matrix[gene_ID,],gene_cluster) # merge together the information from original and cluster method

head(trans_cts_cluster)

Agglometarive -AGNES

Alternatively, we can use the agnes function. These functions behave very similarly; however, with the agnes function you can also get the agglomerative coefficient, which measures the amount of clustering structure found (values closer to 1 suggest strong clustering structure).

library(cluster)
ag.fit <- agnes(log2_col.matrix[gene_ID,])
par(mfrow=c(1,1))
plot(ag.fit,which.plots=2,main='agnes - dendogram',cex=.8,cex.main=2,xlab='')

ag.fit$ac

[1] 0.7517836

di.fit <- diana(log2_col.matrix[gene_ID,])
plot(di.fit,which.plots=2,main='diana – dendogram', cex=.8,cex.main=2,xlab='')

par(mfrow=c(1,2))

plot(ag.fit,which.plots=2,main='agnes-dendogram',cex=.8,cex.main=2,xlab='Genes',col="darkblue")
plot(di.fit,which.plots=2,main='diana–dendogram', cex=.8,cex.main=2,xlab='Genes',col="darkblue")

NA
NA

library(cluster)
library(RColorBrewer)

hmcol2 <- colorRampPalette(brewer.pal(10,"Blues"))(256)

Warning: n too large, allowed maximum for palette Blues is 9
Returning the palette you asked for with that many colors

csc <- rep(hmcol2[50],ncol(log2_col.matrix[gene_ID,]))
heatmap(log2_col.matrix[gene_ID,],scale="column",col=hmcol2,ColSideColors=csc,cexCol=1,cexRow=.5)

Heatmap(log2_col.matrix[gene_ID,],row_split = factor(Colon$Y[gene_ID]),cluster_columns = ,border=T,heatmap_legend_param = list(title = ""), cluster_row_slices = F)

Kmeans cluster

##kmeans

colon_kmeans <- kmeans(log2_col.matrix[gene_ID,], centers = 4, nstart = 25)
fviz_cluster(colon_kmeans, data = log2_col.matrix[gene_ID,])

k2 <- kmeans(log2_col.matrix[gene_ID,], centers = 2, nstart = 25)
k3 <- kmeans(log2_col.matrix[gene_ID,], centers = 3, nstart = 25)
k4 <- kmeans(log2_col.matrix[gene_ID,], centers = 4, nstart = 25)
k5 <- kmeans(log2_col.matrix[gene_ID,], centers = 5, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point",  data = log2_col.matrix[gene_ID,]) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = log2_col.matrix[gene_ID,]) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = log2_col.matrix[gene_ID,]) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = log2_col.matrix[gene_ID,]) + ggtitle("k = 5")

library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

The following object is masked from ‘package:Biobase’:

    combine

The following object is masked from ‘package:BiocGenerics’:

    combine

grid.arrange(p1,p2, p3, p4, nrow = 2)

#Determining Optimal Clusters

set.seed(123)

fviz_nbclust(log2_col.matrix[gene_ID,], kmeans, method = "wss")

library(cluster)
library(factoextra)

fviz_nbclust(log2_col.matrix[gene_ID,], kmeans, method = "silhouette")

# compute gap statistic
set.seed(123)
gap_stat <- clusGap(log2_col.matrix[gene_ID,], FUN = kmeans, nstart = 25,K.max = 10, B = 50)

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 50)  [one "." per sample]:
.................................................. 50 

# Print the result
#print(gap_stat, method = "firstmax")

fviz_gap_stat(gap_stat)

p1<-fviz_nbclust(log2_col.matrix[gene_ID,], kmeans, method = "wss")
p2<-fviz_nbclust(log2_col.matrix[gene_ID,], kmeans, method = "silhouette")
p3<-fviz_gap_stat(gap_stat)

grid.arrange(p1,p2, p3, nrow = 1)

trial<-kmeans(log2_col.matrix[gene_ID,],centers=2)
#table(trial$cluster,Colon$gene.names[gene_ID])

library(fpc)
trial2<-kmeansruns(log2_col.matrix[gene_ID,],krange=2:5,criterion="asw")
trial2$bestk

[1] 2

#table(trial2$cluster,Colon$gene.names[gene_ID])

trial_pam<-pam(log2_col.matrix[gene_ID,],k=2)
#table(trial_pam$cluster,Colon$gene.names[gene_ID])

Question 5

Dimensionality reduction. We decided to reduce the number of samples used based on the correlation to Dim 1 as the one which has the highest impact on explaining the variance.

which(res.pca$var$cor[,1]>=0.85)# which samples have correlation in PCA greater or equal to 0.9 for Dimenssion 1

 1  5  7 13 15 17 25 27 31 32 33 34 35 40 41 51 52 53 
 1  5  7 13 15 17 25 27 31 32 33 34 35 40 41 51 52 53 

# 1  5  7 13 15 17 25 27 31 32 33 34 35 40 41 51 52 53 

samp_ID<-c(1 , 5,  7, 13, 15, 17, 25, 27, 31, 32, 33, 34, 35, 40, 41, 51, 52, 53 )
# length(samp_ID)*100/62 is the % of samples used after reduction to create clusters. Almost 30% of samples is used and the performance is higher.
Colon_subset<-log2_col.matrix[gene_ID,]
Dim_red<-Colon_subset[,samp_ID] # subset of features (samples) from 62 now we have 18
head(Dim_red)

           1         5         7       13       15        17        25        27
26 11.640432 11.710162 11.605570 12.09815 12.65107 12.430876 13.139931 11.255110
31 10.682606 11.370833 11.099395 11.13129 11.61965 11.476299 12.284663 10.492534
43 10.911689 11.259503 10.552497 11.03785 11.44814 11.537175 12.016171 10.298512
47 10.733909 10.892490 10.809409 11.37828 11.97257 11.371233 12.040828 10.276976
62  9.941739 11.434114 10.568063 11.36080 11.02489 10.198126 11.682020  9.946427
66  9.931131  9.674795  9.526905 10.24609 10.91056  9.159745  9.511503  9.346569
         31        32        33        34        35        40        41        51
26 12.63412 11.692359 12.521866 12.458045 11.665094 12.662580 11.657195 12.469209
31 11.43871 10.778548 11.758441 11.468805 10.769609 11.781428 10.394798 11.962385
43 11.35733 10.601648 11.444002 11.373157 10.083139 11.601210 10.332042 11.017357
47 11.84717 10.640656 11.318007 11.937338 10.470172 12.124319 10.307027 10.155285
62 11.16823  9.958680 11.002820 11.497630  9.983184 11.101130  9.986109 10.735953
66 10.58244  9.573926  9.425481  9.873661  9.219780  9.922843  8.905541  9.741088
         52       53
26 13.55813 12.56184
31 12.77695 11.52969
43 12.62469 11.53503
47 12.62698 11.39854
62 12.26840 10.91591
66 11.66462 10.32727

boxplot(t(Dim_red),xlab="Genes",ylab="Expression values",main="Dim reduction subset")# boxplot of all genes

my_set = ExpressionSet(assayData=Dim_red,phenodata=Colon$Y[samp_ID])
plotDensities(my_set,legend=FALSE, main="")

library(cluster)
dim(Dim_red)

[1] 143  18

ag.fit <- agnes(Dim_red)
par(mfrow=c(1,1))
plot(ag.fit,which.plots=2,main='agnes - dendogram',cex=.8,cex.main=2,xlab='')

ag.fit$ac

[1] 0.7864747

di.fit <- diana(Dim_red)
plot(di.fit,which.plots=2,main='diana – dendogram', cex=.8,cex.main=2,xlab='')

par(mfrow=c(1,2))

plot(ag.fit,which.plots=2,main='agnes-dendogram',cex=.8,cex.main=2,xlab='Genes',col="darkblue")
plot(di.fit,which.plots=2,main='diana–dendogram', cex=.8,cex.main=2,xlab='Genes',col="darkblue")

k2 <- kmeans(Dim_red, centers = 2, nstart = 25)
k3 <- kmeans(Dim_red, centers = 3, nstart = 25)
k4 <- kmeans(Dim_red, centers = 4, nstart = 25)
k5 <- kmeans(Dim_red, centers = 5, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point",  data = Dim_red) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = Dim_red) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = Dim_red) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = Dim_red) + ggtitle("k = 5")

library(gridExtra)
grid.arrange(p1,p2, p3, p4, nrow = 2)

p1<-fviz_nbclust(Dim_red, kmeans, method = "wss")
p2<-fviz_nbclust(Dim_red, kmeans, method = "silhouette")
set.seed(123)
gap_stat <- clusGap(Dim_red, FUN = kmeans, nstart = 25,K.max = 10, B = 50)

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 50)  [one "." per sample]:
.................................................. 50 

p3<-fviz_gap_stat(gap_stat)

grid.arrange(p1,p2, p3, nrow = 1)

library(ComplexHeatmap)
Heatmap(as.matrix(t(Dim_red)),row_split = factor(Colon$Y[samp_ID]),cluster_columns = ,border=T,heatmap_legend_param = list(title = ""), cluster_row_slices = F)

heatmap(t(Dim_red))

res.pca <- PCA(Dim_red, graph = FALSE)

fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

fviz_pca_ind(res.pca, col.ind = "contrib", pointsize = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", # Variables color
                col.ind = "#696969"  # Individuals color
                )

End!

Eralda Gjika

https://github.com/EGjika

https://al.linkedin.com/in/eralda-dhamo-gjika-71879128